#DESCRIPTION:
#Barton and Bandis Failure Criterion and FS analysis


#1. libraries and modules

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
#from scipy.stats import beta,gamma,gaussian_kde, kstest, ks_2samp, ttest_ind,norm
#from lmfit import minimize, Parameters,report_errors
#from scipy.interpolate import UnivariateSpline



def BBeq(JRC,JCS,sig_n,phi_r,Ln):
    #equation for Barton-Bandis failure criterion
    JRC,JCS=BB_red(JRC,JCS,Ln)
    return sig_n*np.tan(np.radians(JRC*np.log10(JCS/sig_n)+phi_r))

def BB_red(JRC0,JCS0,Ln):
    #equation for scale-based reduction of JRC, JCS
    L0=0.1 #m
    JRCn=JRC0*(Ln/L0)**(-0.02*JRC0)
    JCSn=JCS0*(Ln/L0)**(-0.03*JRC0)
    return JRCn,JCSn

def generate_random_BB_vars(sign_min,sign_max,
                            phir_min,phir_max,
                            JRC0_min, JRC0_max,
                            JCS0_min,JCS0_max,
                            numrand):
    np.random.seed(12345)
    sig_n=np.random.rand(numrand)*(sign_max-sign_min)+sign_min  
    phi_r=np.random.randint(phir_min,phir_max+1,numrand) 
    JRC0=np.random.randint(JRC0_min,JRC0_max+1,numrand)
    JCS0=np.random.randint(JCS0_min,JCS0_max+1,numrand) 
    return sig_n,phi_r,JRC0,JCS0


def normalize(var):
    #general normalization of array of variables from 0 to 1.
    return (var-min(var))/(max(var)-min(var)*1.)

def BB_sensitivity():
    #0.

    #defining random variables

    numrand=500
    sign_min,sign_max=0.,1.     #MPa
    phir_min,phir_max=10,25     #degrees
    JRC0_min,JRC0_max=0,20     
    JCS0_min,JCS0_max=10,80     #MPa

    sign_max_list=[0.1,.5,1,5.]
    iS=0
    for sign_max in sign_max_list:

        sig_n,phi_r,JRC0,JCS0=generate_random_BB_vars(sign_min,sign_max,
                                phir_min,phir_max,
                                JRC0_min, JRC0_max,
                                JCS0_min,JCS0_max,
                                numrand)
        #1
        Ln_list=[.1,.5,1.,5]
        num_Ln=len(Ln_list)
        numvars=4
        fig,ax=plt.subplots(nrows=numvars,ncols=num_Ln,sharex=True, sharey=True, figsize=(4,4))
        fig.suptitle('sign_max='+str(sign_max)+' MPa\n',fontsize='small')
        symsize=4
        for r in range(num_Ln):
            Ln=Ln_list[r]
            tau=BBeq(JRC0,JCS0,sig_n,phi_r,Ln)
            psig_n=normalize(sig_n)
            pphi_r=normalize(phi_r)
            pJRC0=normalize(JRC0)
            pJCS0=normalize(JCS0)
            var1=psig_n
            var2=pphi_r
            var3=pJRC0
            var4=pJCS0

            ax[0,r].plot(var1,tau,'.',markersize=symsize)            
            ax[1,r].plot(var2,tau,'.',markersize=symsize)
            ax[2,r].plot(var3,tau,'.',markersize=symsize)
            ax[3,r].plot(var4,tau,'.',markersize=symsize)
            
            if r==0:
                ax[0,r].set_ylabel('tau(sig_n)', fontsize='x-small')
                ax[0,r].set_ylim(0,max(tau))
                ax[1,r].set_ylabel('tau(phi_r)', fontsize='x-small')
                ax[2,r].set_ylabel('tau(JRC0)', fontsize='x-small')
                ax[3,r].set_ylabel('tau(JCS0)', fontsize='x-small')

            

            ax[0,r].set_title('Ln='+str(Ln_list[r])+' m', fontsize='x-small')

            ax[0,r].set_yticklabels([])
            ax[1,r].set_yticklabels([])
            ax[2,r].set_yticklabels([])
            ax[3,r].set_yticklabels([])

            ax[0,r].set_xticklabels([])
            ax[1,r].set_xticklabels([])
            ax[2,r].set_xticklabels([])
            ax[3,r].set_xticks([0,1])
            ax[3,r].set_xticklabels([0,1], fontsize='x-small')
        plt.savefig('/home/eglsais/Dropbox/Thesis/Thesis chapters/lyx figures/BB_senA_'+str(iS)+'.png.')           
        plt.show()
        iS=iS+1
##            
    #2.
##    Ln_list=[0.1,0.2,0.5,1]
##    sign_max_list=[0.5,1.,3.,5.]
##    iL=0
##    for Ln in Ln_list:
##        iS=0
##        for sign_max in sign_max_list:
##            if Ln==0.1 and sign_max==0.5:
##                plot_legend=1
##            else:
##                plot_legend=0
##                
##            sig_n,phi_r,JRC0,JCS0=generate_random_BB_vars(sign_min,sign_max,
##                            phir_min,phir_max,
##                            JRC0_min, JRC0_max,
##                            JCS0_min,JCS0_max,
##                            numrand)
##            numvars=3
##            fig,ax=plt.subplots(nrows=numvars,ncols=1,sharex=True, sharey=True,figsize=(3,4))
##            sym='o'
##            symsize=4
##            #phi_r
##            phi_r_list=[10,15,20,25]
##            num_phi_r=len(phi_r_list)
##            for r in range(num_phi_r):
##                phi_r2=phi_r_list[r]
##                tau=BBeq(JRC0,JCS0,sig_n,phi_r2,Ln)
##
##                psig_n=normalize(sig_n)
##                pJRC0=normalize(JRC0)
##                pJCS0=normalize(JCS0)
##                
##                var1=sig_n
##                ax[0].plot(var1,tau,lw=0,marker=sym,markersize=symsize,label='phi_r='+str(phi_r_list[r]))
##                if plot_legend:ax[0].legend(fontsize='x-small',loc='upper left')
##                ax[0].set_ylim(0,max(tau))
##                ax[0].set_yticks([ax[0].get_yticks()[0],ax[0].get_yticks()[len(ax[0].get_yticks())-1]])
##                ax[0].set_yticklabels(ax[0].get_yticks(),fontsize='small')
##                
##
##            #JRC0
##            JRC0_list=[0,6,12,20]
##            num_JRC0=len(JRC0_list)
##            for r in range(num_JRC0):
##                JRC02=JRC0_list[r]
##                tau=BBeq(JRC02,JCS0,sig_n,phi_r,Ln)
##
##                psig_n=normalize(sig_n)
##                pphi_r=normalize(phi_r)
##                pJCS0=normalize(JCS0)
##                
##                var1=sig_n
##                ax[1].plot(var1,tau,lw=0,marker=sym,markersize=symsize,label='JRC0='+str(JRC0_list[r]))
##                if plot_legend:ax[1].legend(fontsize='x-small',loc='upper left')
##                ax[1].set_ylabel('tau, MPa',fontsize='small')
##                ax[1].set_yticks([ax[1].get_yticks()[0],ax[1].get_yticks()[len(ax[1].get_yticks())-1]])
##                ax[1].set_yticklabels(ax[1].get_yticks(),fontsize='small')
##
##            #JCS0
##            JCS0_list=[10,35,65,90]
##            num_JCS0=len(JCS0_list)
##            for r in range(num_JCS0):
##                JCS02=JCS0_list[r]
##                tau=BBeq(JRC0,JCS02,sig_n,phi_r,Ln)
##
##                psig_n=normalize(sig_n)
##                pphi_r=normalize(phi_r)
##                pJRC0=normalize(JRC0)
##                
##                var1=sig_n
##                ax[2].plot(var1,tau,lw=0,marker=sym,markersize=symsize,label='JCS0='+str(JCS0_list[r]))
##                if plot_legend:ax[2].legend(fontsize='x-small',loc='upper left')
##                ax[2].set_xlabel('sig_n, MPa',fontsize='small')
##
##                print ax[2].get_xticks()[len(ax[2].get_xticks())-1]
##                ax[2].set_xticks([ax[2].get_xticks()[0],ax[2].get_xticks()[len(ax[2].get_xticks())-1]])
##                ax[2].set_xticklabels(ax[2].get_xticks(),fontsize='small')
##                ax[2].set_yticks([ax[2].get_yticks()[0],ax[2].get_yticks()[len(ax[2].get_yticks())-1]])
##                ax[2].set_yticklabels(ax[2].get_yticks(),fontsize='small')
##                
##                
##
##            fig.subplots_adjust(left=0.2,top=0.95,hspace=0.1)
##            fig.suptitle('Ln='+str(Ln)+'m',fontsize='medium')
##            #fig.tight_layout()
##            plt.savefig('/home/eglsais/Dropbox/Thesis/Thesis chapters/lyx figures/BB_senB_'+str(iL)+str(iS)+'.png.')           
##            #plt.show()
##            iS=iS+1
##        iL=iL+1
                        




def wedge(beta,base,alpha):
    #computing wedge geometry of slope for barton-bandis stability analysis

    #checking input parameters    
    try:
        beta.shape
    except:
        beta=np.asarray([beta])
    try:
        alpha.shape
    except:
        alpha=np.asarray([alpha])
    if alpha.shape!=beta.shape:
        print "Error: alpha and beta not the same shape"
        return    
    
    H=base*np.tan(np.radians(beta))
    maxx=base*1.5*np.ones(beta.shape) 
    alpha_x_at_H=H/np.tan(np.radians(alpha))
    A_alpha=np.where(alpha_x_at_H>maxx,
                     0.5*maxx*(maxx*np.tan(np.radians(alpha))),
                     0.5*alpha_x_at_H*H)
    A_beta=np.where(alpha_x_at_H>maxx,
                    (0.5*base*H) + H*0.5*base,
                    (0.5*base*H) + H*(alpha_x_at_H-base))
    l_alpha=np.where(alpha_x_at_H>maxx,
                     maxx/np.cos(np.radians(alpha)),
                     alpha_x_at_H/np.cos(np.radians(alpha)))
    return A_beta,A_alpha,l_alpha,alpha,beta

def BB_stability_analysis(k_uw,k_base,k_Ln,
                          v_beta,v_alpha,
                          r_phi_r,r_JRC0,r_JCS0,plot_sensitivity):
    #INPUTS
    #k_* type: float
    #v_* type: float, or if array, must have the same size
    #r_  type: array of random variables having the same size

    #RETURNS
    #df dataframe object containing results of stability analysis, each row assigned to a unique alpha-beta pair

    #computing wedge geometry
    try:
        A1,A2,disc_len,alpha,beta=wedge(v_beta,k_base,v_alpha)
    except:
        return pd.DataFrame()
    
    #creating dataframe
    df=pd.DataFrame()
    df['d_alpha']=alpha
    df['d_beta']=beta
    df['w_Area']=A1-A2
    df['w_Baselength']=disc_len
    df['w_Weight']=df.w_Area*k_uw
    df['w_Sign']=df.w_Weight*np.cos(np.radians(df.d_alpha))/df.w_Baselength
    df['w_Tau']=df.w_Weight*np.sin(np.radians(df.d_alpha))/df.w_Baselength

    #FS analysis
    FSmean_list=np.empty(beta.shape)
    FSstd_list=np.empty(beta.shape)
    FSpcF_list=np.empty(beta.shape)

    for i in range(len(df)):
        tau=BBeq(r_JRC0,r_JCS0,df.w_Sign.values[i],r_phi_r,k_Ln)
        phi=np.degrees(np.arctan(tau/df.w_Sign.values[i]))
        FS=df.w_Sign.values[i]*np.tan(np.radians(phi))/df.w_Tau.values[i]
        FSmean_list[i]=np.mean(FS)
        FSstd_list[i]=np.std(FS)
        FSpcF_list[i]=len(np.where(FS<1)[0])*100./len(FS)

        if plot_sensitivity:
            BB_sensitivity_FS(r_JRC0,r_JCS0,r_phi_r,FS,df.d_beta.values[i],df.d_alpha.values[i])

        

    #appending FS results to dataframe
    df['FS_mean']=FSmean_list
    df['FS_std']=FSstd_list
    df['FS_pcF']=FSpcF_list

    return df        

def BB_sensitivity_FS(r_JRC0,r_JCS0,r_phi_r,FS,beta,alpha):
    fig,ax=plt.subplots(nrows=3,ncols=1,sharey=True)
    nJRC=normalize(r_JRC0)
    nJCS=normalize(r_JCS0)
    nphir=normalize(r_phi_r)

    ax[0].plot(nphir,FS,'.')
    ax[0].set_ylabel('FS(phi_r)')
    ax[1].plot(nJRC,FS,'.')
    ax[1].set_ylabel('FS(JRC0)')
    ax[2].plot(nJCS,FS,'.')
    ax[2].set_ylabel('FS(JCS0)')

    fig.suptitle('beta='+str(beta)+' alpha='+str(alpha))
    fig.tight_layout()
    plt.show()


def sample_analysis():
    #SAMPLE STABILITY ANALYSIS

    #defining constants
    uw=24*10**-3 #MN/m3
    base=30
    Ln=5

    #defining slope gradient (beta) and discontinuity dip (alpha)
    beta=np.asarray([75,75,75,75,75,
                     60,60,60,60,
                     45,45,45,
                     30]) # or import from csv
    alpha=np.asarray([60,50,40,30,20,
                      50,40,30,20,
                      40,30,20,
                      20]) # or import from csv


    #defining random variables

    numrand=1000
    sign_min,sign_max=0.,3.     #MPa
    phir_min,phir_max=10,25     #degrees
    JRC0_min,JRC0_max=0,20     
    JCS0_min,JCS0_max=10,60     #MPa

    sig_n,phi_r,JRC0,JCS0=generate_random_BB_vars(sign_min,sign_max,
                                                  phir_min,phir_max,
                                                  JRC0_min, JRC0_max,
                                                  JCS0_min,JCS0_max,
                                                  numrand)

    df=BB_stability_analysis(uw,base,Ln,
                          beta,alpha,
                          phi_r,JRC0,JCS0,0)

    print df


BB_sensitivity()

    



